In [ ]:
# I load the needed libraries
library(dplyr)
library(scales)
library(GoFKernel)

library(mvtnorm)
library(gplots)

options(warn=-1)
Attaching package: ‘dplyr’


The following objects are masked from ‘package:stats’:

    filter, lag


The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union


The following objects are masked from ‘package:stats’:

    filter, lag


The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union


Loading required package: KernSmooth

KernSmooth 2.23 loaded
Copyright M. P. Wand 1997-2009


Attaching package: ‘gplots’


The following object is masked from ‘package:stats’:

    lowess


Preparation of the simulation¶

I load the functions from the class file:

In [ ]:
source("class_MCMC.R")

I define then the function that I want to use as output of the MCMCs:

In [ ]:
# Function to sampled from: n-dim gaussian with chosen sigmas and centers
# posterior_g_inhom = function (theta) {

#     sigmas = c(1:length(theta))
#     centers = c(seq(length(theta), 1))

#     product = 1
#     for (i in 1:length(theta)) {
#         product = product * exp(-(theta[i] - centers[i])**2/sigmas[i]**2)
#     }

#     return (product)

# }

cauchy1_gauss2 = function (theta) {

    sigmas = c(2.5, 4.3)
    centers = c(0.4, 9)

    product = 1
    for (i in 1:2) {
        product = product * exp(-(theta[i] - centers[i])**2/sigmas[i]**2)
    }

    product = product * (dcauchy(theta[3], -10, 2) + 4*dcauchy(theta[3], 10, 4))

    return (product)

} 

chosen_function = cauchy1_gauss2

Then I only have to determine the parameters for the initialization = the "hyperparameters" of the simulations

In [ ]:
# The initial parameters are:
init = c(2, 4, 10)
std = diag(1, 3)

N = as.integer(1e5)
burn_in = as.integer(1e4)

print_step = as.integer(1e2)
# print_init = as.integer(1e3)

N_tot = N + burn_in

# For Haario:
epsilon = 0.001

Simulations¶

In [ ]:
# MVTNORM 

# Evaluate then the MCMC
mcmc_g = random_steps_mvtnorm (func_wanted = chosen_function, theta_init = init, n_samples = N_tot, sigma = std, print_accept=TRUE)

# Selecting the sequence after the burn-in
mcmc_g = mcmc_g[burn_in:N, ]

# Plotting the results
show_results(mcmc = mcmc_g, init = init, std = std, step = print_step)
Acceptance rate =  76.43091 %
In [ ]:
# MVTNORM GIBBS

mcmc_g = random_steps_mvtnorm_gibbs (func_wanted = chosen_function, theta_init = init, n_samples = N_tot, sigma = std, print_accept=TRUE)

mcmc_g = mcmc_g[burn_in:N, ]

show_results(mcmc = mcmc_g, init = init, std = std, step = print_step)
Acceptance rate =  83.94091 %
In [ ]:
# # SIMPLE ADAPTIVE

# mcmc_g = random_steps_simple (func_wanted = chosen_function, theta_init = init, n_samples = N_tot, sigma = std, print_accept=TRUE, t_0 = burn_in,
#                                 gamma_function = gamma_series_exp, halved_step = burn_in)

# mcmc_g = mcmc_g[burn_in:N, ]

# show_results(mcmc = mcmc_g, init = init, std = std, step = print_step)
In [ ]:
# # SIMPLE ADAPTIVE GIBBS

# mcmc_g = random_steps_simple_gibbs (func_wanted = chosen_function, theta_init = init, n_samples = N_tot, sigma = std, print_accept=TRUE, t_0 = burn_in,
#                                 gamma_function = gamma_series_exp, halved_step = burn_in)

# mcmc_g = mcmc_g[burn_in:N, ]

# show_results(mcmc = mcmc_g, init = init, std = std, step = print_step)
In [ ]:
# HAARIO

mcmc_g = random_steps_haario (func_wanted = chosen_function, theta_init = init, n_samples = N_tot,
                                sigma = std, print_accept=TRUE, t_0 = burn_in, eps = epsilon)

mcmc_g = mcmc_g[burn_in:N, ]

show_results(mcmc = mcmc_g, init = init, std = std, step = print_step)
Acceptance rate =  16.50455 %
Final mean =  0.4193031 9.023863 5.946447 
Final covariance matrix = 
          [,1]      [,2]       [,3]
[1,]  6.125695  11.89185   7.440331
[2,] 11.891855 295.06024 162.545150
[3,]  7.440331 162.54515 894.842882
In [ ]:
# HAARIO GIBBS

mcmc_g = random_steps_haario_gibbs (func_wanted = chosen_function, theta_init = init, n_samples = N_tot,
                                    sigma = std, print_accept=TRUE, t_0 = burn_in, eps = epsilon)

mcmc_g = mcmc_g[burn_in:N, ]

show_results(mcmc = mcmc_g, init = init, std = std, step = print_step)
Acceptance rate =  34.34333 %
Final mean =  0.3690438 8.992309 6.83135 
Final covariance matrix = 
          [,1]      [,2]       [,3]
[1,]  6.470395  11.43808   8.612761
[2,] 11.438083 293.97645 230.063636
[3,]  8.612761 230.06364 830.000087
In [ ]:
# RAO

mcmc_g = random_steps_AM_rao (func_wanted = chosen_function, theta_init = init, n_samples = N_tot, sigma = std, print_accept=TRUE, t_0 = burn_in,
                                gamma_function = gamma_series_exp, halved_step = burn_in/2)

mcmc_g = mcmc_g[burn_in:N, ]

show_results(mcmc = mcmc_g, init = init, std = std, step = print_step)
Acceptance rate =  43.51091 %
Final mean =  0.3750626 9.02413 6.453105 
Final covariance matrix = 
            [,1]        [,2]       [,3]
[1,]  3.13846802 -0.04642791   0.222659
[2,] -0.04642791  9.32302506  -1.267430
[3,]  0.22265900 -1.26743030 405.974801
In [ ]:
# RAO GIBBS

mcmc_g = random_steps_AM_rao_gibbs (func_wanted = chosen_function, theta_init = init, n_samples = N_tot, sigma = std, print_accept=TRUE, t_0 = burn_in,
                                gamma_function = gamma_series_exp, halved_step = burn_in/2)

mcmc_g = mcmc_g[burn_in:N, ]

show_results(mcmc = mcmc_g, init = init, std = std, step = print_step)
Acceptance rate =  59.94212 %
Final mean =  0.4077146 9.051977 6.745132 
Final covariance matrix = 
             [,1]        [,2]        [,3]
[1,]  3.690542211 0.002453223  -0.4264274
[2,]  0.002453223 8.570809601   0.5417770
[3,] -0.426427380 0.541777006 269.2037717
In [ ]:
# GLOBAL

mcmc_g = random_steps_global (func_wanted = chosen_function, theta_init = init, n_samples = N_tot, sigma = std, print_accept=TRUE, t_0 = burn_in,
                                gamma_function = gamma_series_exp, halved_step = burn_in)

mcmc_g = mcmc_g[burn_in:N, ]

show_results(mcmc = mcmc_g, init = init, std = std, step = print_step)
Acceptance rate =  55.27545 %
Final mean =  0.556769 9.134921 5.749301 
Final lambda =  -0.5730512 
Final covariance matrix = 
           [,1]       [,2]       [,3]
[1,]  2.8205581 -0.2656854  -1.682357
[2,] -0.2656854  8.9788965   1.576610
[3,] -1.6823573  1.5766099 219.212919
In [ ]:
# GLOBAL GIBBS

mcmc_g = random_steps_global_gibbs (func_wanted = chosen_function, theta_init = init, n_samples = N_tot, sigma = std, print_accept=TRUE, t_0 = burn_in,
                                gamma_function = gamma_series_exp, halved_step = burn_in)

mcmc_g = mcmc_g[burn_in:N, ]

show_results(mcmc = mcmc_g, init = init, std = std, step = print_step)
Acceptance rate =  69.36121 %
Final mean =  0.1872337 9.060671 6.700329 
Final lambda =  -0.9273487 
Final covariance matrix = 
            [,1]       [,2]         [,3]
[1,]  3.52341291 -0.1182699  -0.02448046
[2,] -0.11826987  8.3311801   0.98728850
[3,] -0.02448046  0.9872885 138.62285448